A Gpu-based Implementation of a Cone Convex Complementarity Approach for Simulating Rigid Body Dynamics with Frictional Contact

In the context of simulating the frictional contact dynamics of large systems of rigid bodies, this paper reviews a novel method for solving large cone complementarity problems by means of a fixed-point iteration algorithm. The method is an extension of the Gauss-Seidel and Gauss-Jacobi methods with overrelaxation for symmetric convex linear complementarity problems. Convergent under fairly standard assumptions, the method is implemented in a parallel framework by using a single instruction multiple data computation paradigm promoted by the Compute Unified Device Architecture library for graphical processing unit programming. The framework supports the analysis of problems with a large number of rigid bodies in contact. Simulation thus becomes a viable tool for investigating the dynamics of complex systems such as ground vehicles running on sand, powder composites, and granular material flow. Introduction Approximating the time evolution of a multibody system in the presence of friction and contact/impact phenomena through numerical simulation continues to be a challenging task. For instance, results reported in [1] indicate that the most widely used ∗Address all correspondence to this author. commercial software package for multibody dynamics simulation has significant difficulties in handling the simple problem of a collection of balls falling in a box, whenever the number of balls becomes larger than 50; in fact, the problem becomes practically intractable when the number of bodies becomes larger than 100. Presented here is an algorithm that can robustly and efficiently approximate the dynamics of rigid bodies undergoing frictional contact [2]. Posing challenges of its own, the case of deformable frictional contact is extensively discussed in [3, 4] and falls outside the scope of this work. Two approaches are most often considered when simulating the dynamics of a multibody system with frictional contact. First is the class of so-called penalty methods, where it is assumed that every time two rigid bodies come in frictional contact, the interaction can be represented by a collection of stiff springs combined with damping elements that act at the interface of the two bodies [5–8]. Implementing these regularization approaches requires little effort beyond that usually associated with developing a multibody dynamics simulation code. Furthermore, this methodology can easily accommodate complex frictional contact mechanisms, as it allows for a large number of “tuning” parameters that, in general, can be adjusted to control the dynamics of the frictional contact interaction. What has prevented the widespread use of this solution is the small step-size at which the 1 Copyright c © 2008 by ASME numerical integration formula, because of stability limitations, is able to advance the simulation. This drawback is related to the stiff spring elements artificially included in the model. Most of the time, this step-size limitation is counterbalanced by the use of implicit integration formulas, a proposition that typically comes at a price as it requires the solution of a discretization nonlinear system at each integration time-step. This in turn leads to a heavy computational burden for scenarios with a large number of active frictional contact events. A second approach, and the one pursued in this work, relies on a different mathematical framework capable of handling applications with hundreds of thousands of frictional contact events. The algorithms in this class draw on time-stepping procedures that produce weak solutions of the differential variational inequality (DVI) that describes the time evolution of rigid bodies with collision, contact, and friction. The DVI as a problem formulation was recently introduced in full generality and classified by differential index [9], though earlier numerical approaches based on DVI formulations do exist [10–12]. Recent work on time-stepping schemes has included both acceleration-force linear complementarity problem (LCP) approaches [13–15] and velocity-impulse LCP-based time-stepping methods [16–19]. The LCPs, obtained as a result of the introduction of inequalities in time-stepping schemes for DVI, coupled with a polyhedral approximation of the friction cone, must be solved at each time step in order to determine the system state configuration as well as the Lagrange multipliers associated with the frictional contact problem [11, 16]. If the simulation entails a large number of contacts and rigid bodies, as is the case of part feeders, packaging machines, and granular flows, the computational burden of classical LCP solvers can become significant. Indeed, a well-known class of numerical solutions for LCPs is based on simplex methods, also known as direct or pivoting methods [20]; however, these methods may exhibit exponential worstcase complexity [21]. They may be impractical even for problems involving as little as a few hundred bodies when friction is present [22, 23]. Further complicating the numerical solution, since the three-dimensional Coulomb friction case leads to a nonlinear complementarity problem (NCP), the use of a polyhedral approximation to morph the NCP into an LCP introduces artificial anisotropy, which affects friction because friction cones become faceted friction pyramids [15–17]. This discrete and finite approximation of friction cones is one of the reasons for the large dimension of the problem that needs to be solved in multibody dynamics with frictional contact. In order to circumvent the limitations imposed by the use of classical LCP solvers and the limited accuracy associated with polyhedral approximations of the friction cone, a parallel fixed-point iteration method with projection on a convex set is proposed, which can directly solve large cone complementarity problems with low computational overhead. The method is based on a time-stepping formulation that solves at every step a cone constrained optimization problem [24]. The time-stepping scheme has been proved to converge in a measure differential inclusion sense to the solution of the original continuous-time DVI. For the proposed approach, about 80% of the computational effort in simulating frictional contact dynamics is spent solving the cone complementarity problem (CCP). The goal of this work is to solve the CCP in parallel by using commodity high-performance computing hardware. Specifically, a methodology is proposed that hinges on the use of parallel computational resources available on NVIDIA’s graphical processing unit (GPU) cards, which can currently handle 12,228 live computational threads simultaneously on the GeForce 8800 series. Tapping into this massively parallel computational resource has been facilitated by NVIDIA’s sharing of a well-integrated application programming interface supported by the Compute Unified Device Architecture (CUDA) library [25]. Formulation of the Multibody Dynamics with Frictional Contact Problem The equations that govern the time evolution of a multibody system can be expressed in the form (see, for instance, [26]) q̇ = L(q)v Mv̇ = fA (t,q,v) , (1) where q = [ r1 ,!1 , . . . ,rnb ,! T nb ]T ∈ R6nb are generalized positions, v = [ ṙ1 , "̄1 , . . . , ṙnb , "̄ T nb ]T ∈ R6nb are generalized velocities, and nb represents the number of bodies in the system. The matrix M is the generalized mass matrix, and fA (t,q,v) represents the vector of generalized applied forces. The convention used here is that any symbol in bold represents a vector or matrix quantity, and an overbar represents a vector quantity represented in the local, body-fixed reference frame associated with a body that is inferred from the context. The formulation of the equations of motion draws on the so-called absolute, or Cartesian, representation of the attitude of each rigid body in the system. For each body j, its orientation is described by a set of three Euler angles, ! j ∈ R3, following the 3-1-3 local rotation sequence (see, for instance, [26]). The rate at which each body changes its orientation is captured by the local angular velocity "̄ j ∈ R3. The location of each body is uniquely determined by a position vector r j = [x j,y j,z j] that specifies where the body-fixed centroidal reference frame is located. The translational velocity of the body is simply ṙ j, where an overdot represents time differentiation. With this set of generalized coordinates, the mass matrix M remains constant and diagonal between any realigning of a body-fixed centroidal reference frame, which can potentially be employed to avoid Euler angles singularities. Also note that, since for each body j there is a locally nonsingular matrix B(! j) such that "̄ j = B(! j)!̇ j , the 2 Copyright c © 2008 by ASME operator L(q) that relates the time derivative of the level-zero generalized coordinates to the level-one generalized coordinates is generally not the identity matrix. Note that no bilateral constraints are present in the current formulation. This case is discussed in [2, 27], and a paper presenting a parallel methodology for the general case of bilateral and unilateral constraints is forthcoming. Two rigid bodies should not penetrate, and, if they are in contact, there should be friction acting at the interface. In order to enforce the nonpenetration constraint, a gap function #(q) ∈ R is assumed to exist and satisfy